uav_lidar_inventory pt. 3

Author

Teja Kattenborn (teja.kattenborn@geosense.uni-freiburg.de); Maximilian Fabi (maximilian.fabi@geosense.uni-freiburg.de)

Published

December 17, 2024

LiDAR Forest Inventory - Part 3

Main contents of this exercise

  • Compare Terrestrial Laser Scanning (TLS)-based and Unmanned Aerial Vehicle (UAV)-based point clouds for forest inventries.

  • Compare UAV and TLS-based instance segmentations to a full forest inventory of the ECOSENSE site (Ettenheim, Germany).

  • Instance segmentations delineate individual trees in point clouds. Here, we test SegmentAnyTree. A small visual recap on SegmentAnyTree, a tool to segment individual trees in LiDAR point clouds across different LiDAR data types:

  • Moreover, we will look into AI-based tree species recognition that was applied on the invidiual tree segmentations.

Remember that we looked into how the point clouds of ULS and TLS differ from one another, here is a quick reminder which shows the two side by side:

Load Packages, paths, etc…

#install.packages(c("lidR","terra", "dplyr", "viridis))
library(lidR)
library(viridis)
library(dplyr)
library(terra)

Reminder: this script assumes that your Quarto document is in the same directory as your data. If this is not the case, consider to set your working directory with setwd("path/to/your/data/folder/").

Load & inspect the data

Inventory data

Load the tree inventory data from ECOSENSE (full inventory).
*.gpkg stands for geopackage format. We can load that in R using terra as SpatVector object.

tree_inv <- vect("tree_inventory_final2.gpkg")
head(tree_inv)
  species tls_name      tls_TH     tls_DBH       tls_X       tls_Y
1      BE    T2792 27.28962149 0.088146948  416700.687 5346611.676
2      BE    T2756 27.65825122  0.16474044  416698.162 5346611.911
3      BE    T2763 26.53824404 0.136313621 416699.4069 5346613.154
4      BE    T2793 28.58093068 0.363458403 416704.5793 5346614.498
5      BE    T2816 28.73643552 0.387836247 416703.9667 5346609.953
6      BE                 <NA>        <NA>        <NA>        <NA>
      tls_Longitude      tls_Latitude odk_add_tree.loop_time
1 7.877502550380541 48.26723262917562   2024-11-26T11:23:27Z
2 7.877468485304306 48.26723441090238   2024-11-26T11:36:16Z
3 7.877485012699943 48.26724575547465   2024-11-26T11:41:34Z
4 7.877554434593113 48.26725852509873   2024-11-26T11:43:50Z
5 7.877547076487936 48.26721756208074   2024-11-26T11:46:38Z
6              <NA>              <NA>   2024-11-26T11:48:11Z
  odk_add_tree.plot_id odk_add_tree.tree_id odk_add_tree.new_tree
1                 <NA>                 <NA>                  true
2                   16                 <NA>                  true
3                 <NA>                 <NA>                  true
4                 <NA>                 <NA>                  true
5                 <NA>                 <NA>                  true
6                 <NA>                 <NA>                  true
                           odk_qrcode.qrcode_id odk_qrcode.qrcode_bu_id
1 https://dt.unr.uni-freiburg.de/ecosense/16_27                        
2                                      05865987                   16_28
3 https://dt.unr.uni-freiburg.de/ecosense/16_29                        
4 https://dt.unr.uni-freiburg.de/ecosense/16_30                        
5 https://dt.unr.uni-freiburg.de/ecosense/16_31                        
6 https://dt.unr.uni-freiburg.de/ecosense/16_32                        
  odk_new_tree_details.circ odk_new_tree_details.other_species
1                        28                               <NA>
2                        52                               <NA>
3                        41                               <NA>
4                       115                               <NA>
5                       125                               <NA>
6                        25                               <NA>
       odk_tree_img odk_gps_tree.Latitude odk_gps_tree.Longitude
1 1732620954736.jpg      48.2672236483333            7.877503075
2 1732621157835.jpg      48.2672254283333       7.87747073833333
3 1732621415109.jpg      48.2672363233333       7.87748745333333
4 1732621558170.jpg      48.2672479383333       7.87755160333333
5 1732621683746.jpg      48.2672075166667       7.87754983833333
6 1732621738341.jpg      48.2672061566667       7.87752802833333
  odk_gps_tree.Altitude odk_gps_tree.Accuracy odk_tree_comment
1               526.659                  0.01                 
2               526.557                  0.01   No barcode tag
3               526.663                  0.01                 
4               526.798                  0.01                 
5               526.669                  0.01                 
6               526.925                  0.01                 
                             odk_PARENT_KEY
1 uuid:e209a154-c01c-4534-b144-cfa7d7713b22
2 uuid:e209a154-c01c-4534-b144-cfa7d7713b22
3 uuid:e209a154-c01c-4534-b144-cfa7d7713b22
4 uuid:e209a154-c01c-4534-b144-cfa7d7713b22
5 uuid:e209a154-c01c-4534-b144-cfa7d7713b22
6 uuid:e209a154-c01c-4534-b144-cfa7d7713b22
                                            odk_KEY       odk_diameter
1 uuid:e209a154-c01c-4534-b144-cfa7d7713b22/loop[1] 0.0891267681314614
2 uuid:e209a154-c01c-4534-b144-cfa7d7713b22/loop[2]  0.165521140815571
3 uuid:e209a154-c01c-4534-b144-cfa7d7713b22/loop[3]  0.130507053335354
4 uuid:e209a154-c01c-4534-b144-cfa7d7713b22/loop[4]  0.366056369111359
5 uuid:e209a154-c01c-4534-b144-cfa7d7713b22/loop[5]  0.397887357729738
6 uuid:e209a154-c01c-4534-b144-cfa7d7713b22/loop[6] 0.0795774715459477
  odk_shift_distance_m odk_shifted_longitude odk_shifted_latitude
1     1.08912676813146              7.877502             48.26723
2     1.16552114081557              7.877468             48.26723
3     1.13050705333535              7.877485             48.26725
4     1.36605636911136              7.877555             48.26726
5     1.39788735772974              7.877547             48.26722
6     1.07957747154595              7.877526             48.26722
        dbh_difference     new_longitude      new_latitude
1 9.798201314614036e-4 7.877502550380541 48.26723262917562
2 7.807008155710227e-4 7.877468485304306 48.26723441090238
3  0.00580656766464599 7.877485012699943 48.26724575547465
4  0.00259796611135904 7.877554434593113 48.26725852509873
5 0.010051110729737966 7.877547076487936 48.26721756208074
6                 <NA>  7.87752802833322  48.2672159464496
                               geometry last_trees TreeID_test_area
1 c(416700.727386307, 5346611.77496821)                       16_27
2 c(416698.331770989, 5346612.08489961)                    05865987
3 c(416699.589345929, 5346613.24243561)                       16_29
4 c(416704.372125069, 5346614.70123457)                       16_30
5 c(416704.175926162, 5346610.24241929)                       16_31
6  c(416702.550551442, 5346609.7940927)                       16_32
                      layer                                   path TreeID
1 last_trees_1 — last_trees last_trees_1.gpkg|layername=last_trees  16_27
2 last_trees_1 — last_trees last_trees_1.gpkg|layername=last_trees  16_28
3 last_trees_1 — last_trees last_trees_1.gpkg|layername=last_trees  16_29
4 last_trees_1 — last_trees last_trees_1.gpkg|layername=last_trees  16_30
5 last_trees_1 — last_trees last_trees_1.gpkg|layername=last_trees  16_31
6 last_trees_1 — last_trees last_trees_1.gpkg|layername=last_trees  16_32
     species_name
1 Fagus_sylvatica
2 Fagus_sylvatica
3 Fagus_sylvatica
4 Fagus_sylvatica
5 Fagus_sylvatica
6 Fagus_sylvatica

….puh, many attributes. For sure not all of them are useful for this analysis. We will mostly focus on the species, the DBH and of course the geolocation:

plot(tree_inv, col = "red")

# please ignore
#tree_inv_spec <- read.csv("lookup_odk2.csv", sep = ",", header = T)
#tree_inv2 <- merge(tree_inv, tree_inv_spec, by.x = "species", by.y = "abbr")
#writeVector(tree_inv2, "D:/output_file.gpkg", filetype = "GPKG")

TLS & UAV-based inventory data

UAV-based point clouds with AI-based tree segmentation

pc_uav <- readLAS("uav_pointcloud_red.las")
pc_uav
class        : LAS (v1.4 format 6)
memory       : 326.8 Mb 
extent       : 416648.3, 416838.9, 5346560, 5346871 (xmin, xmax, ymin, ymax)
coord. ref.  : NA 
area         : 41168 units²
points       : 6.12 million points
density      : 148.66 points/units²

TLS-based point clouds with AI-based tree segmentation

pc_tls <- readLAS("tls_pointcloud_red.las")
pc_tls
class        : LAS (v1.4 format 6)
memory       : 294.8 Mb 
extent       : 416648.2, 416839, 5346560, 5346871 (xmin, xmax, ymin, ymax)
coord. ref.  : NA 
area         : 41140 units²
points       : 5.52 million points
density      : 134.17 points/units²

Let´s plot the data for a visual comparison:

col <- sample(viridis(2070, option = "mako"))
p <- plot(pc_uav, color = "PredInstance", pal = col, legend = F, axis = T, nbreaks = length(col))
col <- sample(viridis(2070, option = "mako"))
p <- plot(pc_tls, color = "PredInstance", pal = col, legend = F, axis = T, nbreaks = length(col))

UAV-based segmentation vs full inventory

Remember how we extracted the tree cordinates in the previous exercise? We will apply this again.

Normalize point cloud on Z-axis

mycsf <- csf(TRUE, 1, 1, time_step = 1)
pc_uav <- classify_ground(pc_uav, mycsf)
pc_uav_n <- normalize_height(pc_uav, knnidw())
Inverse distance weighting: [=====================================-------------] 74% (8 threads)
Inverse distance weighting: [=====================================-------------] 75% (8 threads)
Inverse distance weighting: [======================================------------] 76% (8 threads)
Inverse distance weighting: [======================================------------] 77% (8 threads)
Inverse distance weighting: [=======================================-----------] 78% (8 threads)
Inverse distance weighting: [=======================================-----------] 79% (8 threads)
Inverse distance weighting: [========================================----------] 80% (8 threads)
Inverse distance weighting: [========================================----------] 81% (8 threads)
Inverse distance weighting: [=========================================---------] 82% (8 threads)
Inverse distance weighting: [=========================================---------] 83% (8 threads)
Inverse distance weighting: [==========================================--------] 84% (8 threads)
Inverse distance weighting: [==========================================--------] 85% (8 threads)
Inverse distance weighting: [===========================================-------] 86% (8 threads)
Inverse distance weighting: [===========================================-------] 87% (8 threads)
Inverse distance weighting: [============================================------] 88% (8 threads)
Inverse distance weighting: [============================================------] 89% (8 threads)
Inverse distance weighting: [=============================================-----] 90% (8 threads)
Inverse distance weighting: [=============================================-----] 91% (8 threads)
Inverse distance weighting: [==============================================----] 92% (8 threads)
Inverse distance weighting: [==============================================----] 93% (8 threads)
Inverse distance weighting: [===============================================---] 94% (8 threads)
Inverse distance weighting: [===============================================---] 95% (8 threads)
Inverse distance weighting: [================================================--] 96% (8 threads)
Inverse distance weighting: [================================================--] 97% (8 threads)
Inverse distance weighting: [=================================================-] 98% (8 threads)
Inverse distance weighting: [=================================================-] 99% (8 threads)
Inverse distance weighting: [==================================================] 100% (8 threads)
plot(pc_uav_n, axis = T)

Cut a horizonal transect for stem extraction (we slice the point cloud close to the ground).

pc_uav_n_trans <- filter_poi(pc_uav_n, Z >= 0.1 & Z <= 4)
col <- sample(viridis(2070, option = "mako"))
plot(pc_uav_n_trans, color = "PredInstance", pal = col, legend = F, axis = T, nbreaks = length(col))

Derive the centroids for stems (as proxy for the tree position):

# Calculate centroids for each instance
centroids_uav <- pc_uav_n_trans@data %>%
  group_by(PredInstance) %>%
  summarize(
    X = mean(X),
    Y = mean(Y),
    Z = 0,
    species = as.integer(names(which.max(table(species_id)))),
    .groups = "drop"  # Avoid unnecessary grouping in the result
  )

# Print centroids
print(centroids_uav)
# A tibble: 1,751 × 5
   PredInstance       X        Y     Z species
          <int>   <dbl>    <dbl> <dbl>   <int>
 1            0 416745. 5346740.     0       0
 2            1 416779. 5346568.     0      19
 3            2 416772. 5346570.     0      10
 4            5 416784. 5346569.     0      10
 5            6 416784. 5346571.     0      10
 6            7 416780. 5346569.     0      10
 7            8 416787. 5346566.     0      25
 8            9 416796. 5346561.     0      19
 9           10 416790. 5346563.     0      10
10           12 416797. 5346568.     0      10
# ℹ 1,741 more rows

Convert the centroids to a geospatial vector format (SpatVector), so we can compare it with the inventory data (which is also in vector format):

centroids_uav_sv <- vect(centroids_uav, geom = c("X", "Y"), crs = "EPSG:32632")

Visually compare the UAV-based centroids with the inventory data

plot(centroids_uav_sv)
points(tree_inv, col = "red")

Quantitatively compare the UAV-based centroids with the inventory data

buffer_tree_inv <- buffer(tree_inv, width = 1)
plot(buffer_tree_inv)
points(tree_inv, col = "red")

matches_uav <- intersect(buffer_tree_inv, centroids_uav_sv)
plot(tree_inv, col = "red")
points(matches_uav, col = "blue", cex = 1.3)

nrow(matches_uav)/nrow(tree_inv)*100
[1] 66.16786

Is there a systematic bias in the detection?

unmatched_uav <- buffer_tree_inv[!(buffer_tree_inv$TreeID %in% matches_uav$TreeID), ]
par(mfrow = c(1,2))
hist(as.numeric(matches_uav$tls_DBH),
     xlim = c(0,0.9),
     ylim=c(0,200),
     main = "DBH matched",
     xlab = "DBH [cm]")
hist(as.numeric(unmatched_uav$tls_DBH),
     xlim = c(0,0.9),
     ylim=c(0,200),
     main = "DBH unmatched",
     xlab = "DBH [cm]")

Not all trees are matched, but at least there seems to be no strong bias for small/large trees.

TLS-based segmentation vs full inventory

Same procedure as for the UAV data. First, we normalize the point cloud…

mycsf <- csf(TRUE, 1, 1, time_step = 1)
pc_tls <- classify_ground(pc_tls, mycsf)
pc_tls_n <- normalize_height(pc_tls, knnidw())
Inverse distance weighting: [=======================================-----------] 79% (8 threads)
Inverse distance weighting: [========================================----------] 80% (8 threads)
Inverse distance weighting: [========================================----------] 81% (8 threads)
Inverse distance weighting: [=========================================---------] 82% (8 threads)
Inverse distance weighting: [=========================================---------] 83% (8 threads)
Inverse distance weighting: [==========================================--------] 84% (8 threads)
Inverse distance weighting: [==========================================--------] 85% (8 threads)
Inverse distance weighting: [===========================================-------] 86% (8 threads)
Inverse distance weighting: [===========================================-------] 87% (8 threads)
Inverse distance weighting: [============================================------] 88% (8 threads)
Inverse distance weighting: [============================================------] 89% (8 threads)
Inverse distance weighting: [=============================================-----] 90% (8 threads)
Inverse distance weighting: [=============================================-----] 91% (8 threads)
Inverse distance weighting: [==============================================----] 92% (8 threads)
Inverse distance weighting: [==============================================----] 93% (8 threads)
Inverse distance weighting: [===============================================---] 94% (8 threads)
Inverse distance weighting: [===============================================---] 95% (8 threads)
Inverse distance weighting: [================================================--] 96% (8 threads)
Inverse distance weighting: [================================================--] 97% (8 threads)
Inverse distance weighting: [=================================================-] 98% (8 threads)
Inverse distance weighting: [=================================================-] 99% (8 threads)
Inverse distance weighting: [==================================================] 100% (8 threads)
plot(pc_tls_n, axis = T)

… create a transect of stems…

pc_tls_n_trans <- filter_poi(pc_tls_n, Z >= 0.1 & Z <= 4)
col <- sample(viridis(2070, option = "mako"))
plot(pc_tls_n_trans, color = "PredInstance", pal = col, legend = F, axis = T, nbreaks = length(col))

… find the centroids of the stems.

# Calculate centroids for each instance
centroids_tls <- pc_tls_n_trans@data %>%
  group_by(PredInstance) %>%
  summarize(
    X = mean(X),
    Y = mean(Y),
    Z = 0,
    species = as.integer(names(which.max(table(species_id)))),
    .groups = "drop"  # Avoid unnecessary grouping in the result
  )

# Print centroids
print(centroids_tls)
# A tibble: 2,093 × 5
   PredInstance       X        Y     Z species
          <int>   <dbl>    <dbl> <dbl>   <int>
 1            0 416740. 5346724.     0       0
 2          797 416718. 5346589.     0      10
 3          798 416723. 5346588.     0      10
 4          802 416727. 5346587.     0      10
 5          803 416735. 5346586.     0      10
 6          807 416736. 5346582.     0      10
 7          808 416730. 5346586.     0      10
 8          809 416737. 5346583.     0       3
 9          810 416734. 5346583.     0      10
10          811 416742. 5346580.     0       3
# ℹ 2,083 more rows

Let´s compare the result to the inventory data:

centroids_tls_sv <- vect(centroids_tls, geom = c("X", "Y"), crs = "EPSG:32632")
plot(centroids_tls_sv)
points(tree_inv, col = "red")

buffer_tree_inv <- buffer(tree_inv, width = 1)
plot(buffer_tree_inv)
points(tree_inv, col = "red")

matches_tls <- intersect(centroids_tls_sv, buffer_tree_inv)
plot(buffer_tree_inv)
points(matches_tls, col = "blue")

nrow(matches_tls)/nrow(tree_inv)*100
[1] 81.52245

TLS & UAV tree height comparison

Above we have looked at the coordinates. But we know that we can extract much more, such as the height.

Extract maximum Z value per instance for both TLS and UAV data:

# Find the highest point for each instance
uav_tree_heights <- pc_uav_n@data %>%
  group_by(PredInstance) %>%
  slice_max(Z, n = 1) %>%  # Get the row with the maximum Z value for each instance
  ungroup()

uav_tree_heights <- uav_tree_heights[, c("X", "Y", "Z", "PredInstance")]

# View the results
print(uav_tree_heights)
# A tibble: 2,161 × 4
         X        Y     Z PredInstance
     <dbl>    <dbl> <dbl>        <int>
 1 416754. 5346752.  46.6            0
 2 416778. 5346567.  28.4            1
 3 416772. 5346569.  24.7            2
 4 416778. 5346569.  28.6            4
 5 416788. 5346570   13.1            5
 6 416788. 5346571.  30.4            6
 7 416780. 5346570.  29.2            7
 8 416788. 5346566.  34.2            8
 9 416796. 5346561.  27.0            9
10 416791. 5346562.  26.7           10
# ℹ 2,151 more rows
# Find the highest point for each instance
tls_tree_heights <- pc_tls_n@data %>%
  group_by(PredInstance) %>%
  slice_max(Z, n = 1) %>%  # Get the row with the maximum Z value for each instance
  ungroup()

tls_tree_heights <- tls_tree_heights[, c("X", "Y", "Z", "PredInstance")]

# View the results
print(tls_tree_heights)
# A tibble: 2,578 × 4
         X        Y     Z PredInstance
     <dbl>    <dbl> <dbl>        <int>
 1 416750. 5346755.  37.9            0
 2 416711. 5346592.  15.7          669
 3 416760. 5346574.  19.0          700
 4 416756. 5346575.  18.3          703
 5 416768. 5346570.  24.4          709
 6 416785. 5346564.  27.6          717
 7 416699. 5346597.  11.8          787
 8 416692. 5346599.  18.6          789
 9 416716. 5346590.  28.7          797
10 416719. 5346589.  26.4          798
# ℹ 2,568 more rows

Let´s compare this visually using histograms:

par(mfrow = c(1,2))
hist(tls_tree_heights$Z, ylim = c(0,900), main = "TLS-based height distribution")
hist(uav_tree_heights$Z, ylim = c(0,900), main = "UAV-based height distribution")

paste("number TLS-based trees: ", nrow(tls_tree_heights))
[1] "number TLS-based trees:  2578"
print("TLS-based tree heights")
[1] "TLS-based tree heights"
summary(tls_tree_heights$Z)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   0.00   10.56   23.44   19.51   27.96   46.18 
paste("number UAV-based trees: ", nrow(uav_tree_heights))
[1] "number UAV-based trees:  2161"
print("UAV-based tree heights")
[1] "UAV-based tree heights"
summary(uav_tree_heights$Z)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  2.065  12.636  25.462  21.584  28.962  48.443 

When looking at the transects … can you explain the results?

AI-based Species identification

Let´s first have look at the reality (forest inventory data of ECOSENSE):

par(mar = c(5, 15, 2, 5))
barplot(sort(table(tree_inv$species_name)), las = 1, horiz = T)

The UAV- and TLS-based predictions come in species-codes, not species names. So we first load a look up table to attach (merge) the actual names to our data:

species_ai <- read.csv("lookup_ai.csv")

… attach the species name to the point clouds….

pc_tls@data <- merge(pc_tls@data, species_ai, by = "species_id", all.x = TRUE)
pc_uav@data <- merge(pc_uav@data, species_ai, by = "species_id", all.x = TRUE)

… and to the extract stem coordinates (centroids):

centroids_tls_sv2 <- merge(centroids_tls_sv , species_ai, by.x ="species", by.y = "species_id", all.x = TRUE)
centroids_uav_sv2 <- merge(centroids_uav_sv , species_ai, by.x ="species", by.y = "species_id", all.x = TRUE)

Visual comparison of the species classification

p <- plot(pc_tls, color = "species_id", legend = T, axis = T)
p <- plot(pc_uav, color = "species_id", legend = T, axis = T)

Okay, the results “look” already quite different. Let´s check the stats:

uav_n_spec <- pc_uav@data %>%
  as.data.frame() %>%
  group_by(species) %>%
  summarise(number_of_instances = n_distinct(PredInstance)) %>%
  arrange(desc(number_of_instances))

uav_n_spec
# A tibble: 15 × 2
   species               number_of_instances
   <chr>                               <int>
 1 Fagus_sylvatica                      1110
 2 Betula_pendula                        408
 3 Pinus_radiata                         185
 4 Carpinus_betulus                       91
 5 Pseudotsuga_menziesii                  83
 6 Quercus_rubra                          77
 7 Pinus_sylvestris                       68
 8 Picea_abies                            63
 9 Acer_pseudoplatanus                    34
10 Acer_saccharum                         16
11 Abies_alba                             12
12 Acer_campestre                          6
13 Quercus_petraea                         4
14 Larix_decidua                           3
15 Quercus_robur                           1
tls_n_spec <- pc_tls@data %>%
  as.data.frame() %>%
  group_by(species) %>%
  summarise(number_of_instances = n_distinct(PredInstance)) %>%
  arrange(desc(number_of_instances))

tls_n_spec
# A tibble: 23 × 2
   species               number_of_instances
   <chr>                               <int>
 1 Fagus_sylvatica                      1024
 2 Acer_saccharum                        626
 3 Pseudotsuga_menziesii                 196
 4 Abies_alba                             95
 5 Picea_abies                            92
 6 Pinus_sylvestris                       56
 7 Quercus_robur                          52
 8 Corylus_avellana                       46
 9 Betula_pendula                         45
10 Carpinus_betulus                       43
# ℹ 13 more rows

sort(table(pc_tls$species))

par(mar = c(10, 5, 2, 5))

# Convert species to a named vector
species_counts <- setNames(tls_n_spec$number_of_instances, tls_n_spec$species)

# Create the barplot
barplot(species_counts,
        las = 2,              # Rotate axis labels for readability
        col = "skyblue",      # Add some color
        main = "Species numbers from TLS data",  # Add a title
        ylab = "Number of Instances")  # Add Y-axis label

par(mar = c(10, 5, 2, 5))

# Convert species to a named vector
species_counts <- setNames(uav_n_spec$number_of_instances, uav_n_spec$species)

# Create the barplot
barplot(species_counts,
        las = 2,              # Rotate axis labels for readability
        col = "skyblue",      # Add some color
        main = "Species numbers from UAV data",  # Add a title
        ylab = "Number of Instances")  # Add Y-axis label

Compare the tree species distribution among datasets:

# Extract unique species across all SpatVectors
all_species <- unique(c(tree_inv$species_name, centroids_tls_sv2$species.y, centroids_uav_sv2$species.y))

tree_inv$species.y <- tree_inv$species_name

# Create a color palette for all species
species_colors <- setNames(rainbow(length(all_species)), all_species)

# Function to assign colors based on species
get_colors <- function(spatvector, color_mapping) {
  color_mapping[spatvector$species.y]
}

# Adjust layout to include space for the legend
par(mfrow = c(1, 3), mar = c(4, 4, 2, 1))  # 2x2 grid, adjust margins as needed


plot(tree_inv, 
     col = get_colors(tree_inv, species_colors), 
     main = "Inventory",
     legend = FALSE)

plot(centroids_tls_sv2, 
     col = get_colors(centroids_tls_sv2, species_colors), 
     main = "Species by TLS",
     legend = FALSE)

plot(centroids_uav_sv2, 
     col = get_colors(centroids_uav_sv2, species_colors), 
     main = "Species by UAV",
     legend = FALSE)

par(mfrow = c(1, 3), mar = c(1, 4, 1, 1))
plot.new()  # Create an empty plot for the legend
legend("center", legend = names(species_colors), fill = species_colors, title = "Species")

Open Questions & Outlook

  • How could we potentially improve the tree species classification?

  • Can we trust UAV-based or TLS-based inventories?

  • Can we trust field-based inventories in terms of coverage, sampling and measurement uncertainty?

  • What are related advantages or disadvantages?

  • How can we move on with UAV- or TLS-based inventories to estimate basal area, biomass or timber volume?

  • How will forest inventories look like in 2050?